qm2ff module

The qm2ff module mainly implement the QM2FF class which implement an automatic method in order to set up a classical force fields from a QM calculation. The methodology is based on the following references:

  1. Bautista, E. J.; Seminario, J. M. Int. J. Quantum Chem. 2008, 108 (1), 180–188.
  2. Seminario, J. M. Int. J. Quantum Chem. 1996, 60 (7), 1271–1277.

Summary of implementation strategy

This is a short summary of the choice made for the implementation.

About bonds

The core part is the definition of the bonds according to a bond search strategy.

  • hessian: bonds are selected if the bond constant is real and positive.
  • bond_order: bonds are selected if the bond order value is not null.
  • hybrid: union of the hessian and bond_order strategy
  • connectivity: bonds are selected according to the provide connectivity.

At the end of the bonds definition, the bond_list list is filled in order to further build the angle, dihedral and improper lists. The bond_list is filled from either the bond_order or the connectivity. That means that the bonds identified from the force constants computed in the hessian bond search strategy are not considered.

About angles

We consider that a bending angle exist between three atoms i, j and k if atom j is bonded to both atoms i and k.

Next, all force constants are considered and a warning is printed if some eigenvalues present an imaginary part.

At the end, the angle_list and the list angles are consistant.

About dihedrals

Two types of dihedrals are available. Parabolic dihedrals are represented by a quadratic potential energy such as:

\[V_{\phi} = k_{\phi}^{para} (\phi - \phi_o)^2\]

harminonic dihedrals are represented by a sinusiidal potential energy such as:

\[V_{\phi} = k_{\phi}^{harm} \left[1 + \cos(p \phi - \delta) \right]\]

where $p$ is the periodicity and $delta$ is the phase of the potential.

The relation between \(k_{\phi}^{para}\) and \(k_{\phi}^{harm}\) reads

\[k_{\phi}^{harm} = \frac{k_{\phi}^{para}}{p^2}\]

A dihedral angle between atoms i, j, k and l correspond to a torsion aroung the bond between atoms j and k. That dihedral angle exists if:

  • the bonds between i-j, j-k and k-l exist
  • if i, j and k are not aligned
  • if j, k and l are not aligned

Next, all force constants are considered and a warning is printed if some eigenvalues present an imaginary part.

Parabolic dihedrals are computed first and harmonic dihedrals, if required, are computed from the parabolic ones. The periodicity is computed following the algorithm suggested by Nilson et al. from the number of neighbors of atoms j and k.

  • if j and k have got the same number of neighbors

    • number of neighbors is 2, periodicity is 1
    • number of neighbors is 3, periodicity is 2
    • number of neighbors is 4, periodicity is 3

  • if j and k have’nt got the same number of neighbors

    • periodicity is 1 if the number of neighbors is 2
    • periodicity is 6 in other cases

If the number of neighbors is larger than 4, until now, an exception is raises. That case might be encounterd in the case of metallic systems.

At the end, the dihedral_list and the list dihedrals, are consistant. A dihedral angle initially in the dihedral_list is removed if it does not satisfy the above conditions.

About impropers

We consider that an improper angle exist between atoms i, j, k and l if

  • atom i is bonded to atoms j, k and l
  • atoms j, k and l are not aligned

Next, all force constants are considered and a warning is printed if some eigenvalues present an imaginary part.

At the end, the improper_list and the list impropers, are consistant. An improper angle initially in the improper_list is removed if it does not satisfy the above conditions.

The QM2FF class

class mammoth.qm2ff.QM2FF(molecule, connectivity=None, cutoffb=2.5, bond_search='hessian', imag_th=0.001, dihedral_type='harmonic', improper_type='all', cutoff_sp2=10.0)[source]

A class that implements the Seminario method to get a force field

  1. Bautista, E. J.; Seminario, J. M. Int. J. Quantum Chem. 2008, 108 (1), 180–188.
  2. Seminario, J. M. Int. J. Quantum Chem. 1996, 60 (7), 1271–1277.

TODO: update the doc with new references when the methods will be improved

General attributes

molecule

A molecule object which contains at least the hessian matrix.

natom

Number of atom in the molecule.

bond_list

List of tuples of two atom indexes which define a bond according to the bond orders or the connectivity.

Warning: depending on the bond search strategy, the bond_list and bonds are not consistents because the bond_list only contains bond defined by the bond_orders or the connectivity.

bonds

A list of Bond objects which describe a bond with atom indexes, bond length, force constants and atomic species.

angle_list

List of tuples of 3 atom indexes which define an angle.

angles

A list of Angle objects which describe an angle with atom indexes, angle value in degree, force constants and atomic species.

dihedral_list

List of tuples of 4 atom indexes which define a dihedral angle.

dihedrals

A list of Dihedral objects which describe a dihedral angle with atom indexes, dihedral angle value in degree, force constants and atomic species. The periodicity can be computed if neighbors are provided.

improper_list

List of tuples of 4 atom indexes which define an improper torsion.

impropers

A list of Improper objects which describe an improper torsionwith atom indexes, torsion angle value in degree, force constants and atomic species.

connectivity

The connectivity of the molecule as a list of atom index pairs. The connectivity depends on the bond search strategy.

Parameters for the calculation

cutoffb

Cutoff distance in angstrom for in order to build the bond list.

Bond search strategy :

  • hessian (default): bond exists if the force constant is real
  • bond_order: bond exists if it is in the molecule.bond_orders dictionary
  • hybrid: bond exists if the force constant is real or if it is in the molecule.bond_orders dictionary
  • connectivity: consider only bonds defined by the connectivity. In that
    case, the connectivity must be provided.
imag_th

Threshold in order to determine if an imaginary part is zero.

dihedral_type

Dihedral type is harmonic or parabolic depending on the functional form of the potential energy. Default is harmonic.

Note about complex eigenvalues : If the sub hessian matrix of used to compute the force constant of each internal coordinates cannot be diagonnalized in R, numpy return the complex eigenvectors and complex eigenvalues. Thus the test concerning the imaginary part of the eigenvalues is relevant even if the hessian matrix is not explicitely a complex matrix.

Set up QM2FF object to compute a force field

Parameters:
  • molecule (Molecule) – A molecule object with hessian matrix and bond orders
  • cutoffb (float) – cutoff distance for bond search
  • bond_serch (str) – Define how bonds are identified. Authorized values are “hessian” (default), “bond_order” or “hybrid”
  • imag_th (float) – threshold under which imaginary part is zero (default 0.001)
  • connectivity (list) – if bond_search is ‘connectivity’ you have to provide the molecular connectivity as in a PDB file: [(0, 1), (1, 2, 3), (2, 5), …] where the first index is bonded to the following indexes.
  • dihedral_type (str) – ‘harmonic’ or ‘parabolic’
  • improper_type (str) – ‘all’ or ‘sp2’
  • cutoff_sp2 (float) – if improper_type is sp2, cutoff …
get_PDB_connectivity()[source]

Return the connectivity in a PDB format. The connectivity depends on the chosen bond search strategy. The output string is such as :

CONNECT 1 3 CONNECT 2 3 CONNECT 3 2 4 1 CONNECT 4 3 7 5 6 CONNECT 5 4 CONNECT 6 4 CONNECT 7 4
get_angle_list()[source]

Build the list of possible angles according to the bond order list or the connectivity. In consequence the angle list does not depend on the bond search strategy.

Angles are store in the angle_list attribute.

get_angles()[source]

Compute the force constant associated to angles store in the angle_list attribute and sotre the results in a list of Angle objects.

get_bonds()[source]

Build the list of bonds according to the bond search strategy and compute the associated force constants. Bond search strategy are :

  • hessian: bond exists if the force constant is real
  • bond_order: bond exists if it is in the molecule.bond_orders dictionary
  • hybrid: bond exists if the force constant is real or if it is in the molecule.bond_orders dictionary
  • connectivity: only consider bond defined by the connectivity

The distance cutoff between atom is defined by self.cutoffb.

get_dihedral_list()[source]

Build the list of possible dihedral angles according to the bond order list or the connectivity. In consequence the dihedral angle list does not depend on the bond search strategy.

Dihedral are store in the dihedral_list attribute.

get_dihedrals_harmonic()[source]

Compute the force constant associated to dihedral angles considering an harmonic expression of the potential energy along the coordinate. The force constant of the harmonic potential is computed from the force constant of the parabolic potential. Assuming an harmonic potential

V_phi = K [1 + cos(n phi - d)]

The force constant K reads

K = K_p / n^2

where K_p is the parabolic force constant.

All dihedral angles stored in the dihedral_list attribute are considered and the force constants are store in a list of Dihedral objects.

get_dihedrals_parabolic()[source]

Compute the force constant associated to dihedral angles considering a parabolic expression of the potential energy along the coordinate as proposed in the Seminario method (IJQC 30, 1996, 1271-1277).

V_phi = 1 / 2 k_phi (phi - phi_0)^2

All dihedral angles stored in the dihedral_list attribute are considered and the force constants are store in a list of Dihedral objects.

get_forcefield()[source]

Compute all parameters and return a ClassicalFrorceField object

get_improper_list()[source]

Build the list of all possible improper torsions according to the bond order list or the connectivity. In consequence the improper list does not depend on the bond search strategy.

Improper torsions are store in the improper_list attribute.

get_impropers()[source]

Compute the force constant associated to improper torsions store in the improper_list attribute and sotre the results in a list of Improper objects.

set_PDB_connectivity(connectivity)[source]

Convert the connectivity initially in a format such as PDB in a list of atom index pairs and store the resulting list in a set object.

WARNING: atom indexes in PDB usually start at 1. Thus atom indexes are translated before being store in the connectivity.

Example of connectivity in PDB format:

CONNECT 1 3 CONNECT 2 3 CONNECT 3 1 2 4 CONNECT 4 3 5 6 7 CONNECT 5 4 CONNECT 6 4 CONNECT 7 4

A list such as the following is expected:

[[1, 3], [2, 3], [3, 1, 2, 4], [4, 3, 5, 6, 7], [5, 4], [6, 4], [7, 4]]
connectivity

Return the connectivity as a list of atom index pairs. If the connectivity was not defined, the returned connectivity is the one obtained from the chosen bond search strategy. In consequence, the connectivity may differ from the bond orders.